;+
; NAME:
;       IRDC_MORPH
;
; PURPOSE:
;       Perform morphological matching between BGPS images and the
;       GLIMPSE 8um images.  This method uses a simple radiative
;       transfer model to create model 8-um images based on BGPS
;       images and the estimated I_MIR.  By using the model of 8-um
;       foreground fraction from the irdc_dist_model/ work, this
;       routine creates model 8-um images as a function of
;       heliocentric distance.  Comparison of these model images with
;       the GLIMPSE image produces a likelihood as a function of
;       distance for a match.
;
; CATEGORY:
;       distance-omnibus IRDC association routine
;
; CALLING SEQUENCE:
;       prob = IRDC_MORPH(struct [,/MAKE_PS])
;
; INPUTS:
;       STRUCT   -- BGPS structure containing coordinates, Bolocat
;                   fitted ellipse information, and others.
;
; OPTIONAL INPUTS:
;       TOPHAT   -- Diameter of the top-hat restriction on the
;                   morphological matching region [arcsec].
;                   Default = 0 (i.e. no top-hat imposed)
;       Td       -- Dust temperature to use for conversion of BGPS
;                   flux density to 8-um optical depth.
;                   [Default = !IRDC.Td]
;       FMEAS    -- Calculated f_meas for object in STRUCT (used for
;                   plotting).
;       MDIST    -- Model fit distance for object in STRUCT (used for
;                   plotting).
;       SUFF     -- Suffix (like "_q8") of the UBC_PROC images to be
;                   read in.  [Default = '']
;
; KEYWORD PARAMETERS:
;       MAKE_PS  -- Create postscript files containing the various
;                   images and plot of likelihood as a function of
;                   distance.  Creates files in the './local/output'
;                   directory.
;       VERBOSE  -- Be verbose with output, including object # and
;                   statistics related to the TAUMAP.
;       CLEM     -- Use the Clemens (1985) rotation curve instead of
;                   the standard Reid (2009) curve.
;       BEAMMASK -- Use a Gaussian (33" FWHM) weighting on the
;                   chi^2 comparison between the model image and the
;                   smoothed GLIMPSE image.
;
; OUTPUTS:
;       PROB     -- Likehihood (or chisqr) as a function of distance
;                   for association between the BGPS image and the
;                   GLIMPSE 8-um image.
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       POSTAGE_BLOCK -- Block holding the postage-stamp images loaded
;                        into memory by load_common_postage.pro.
;                        NOTE: Use of this block is optional, though
;                        it speed up execution.
;       CUBE_BLOCK    -- Block holding the FFORE data cube (l,b,d)
;                        needed to generate the DPDF.  If
;                        omni_load_ffore.pro is not called previous to
;                        this routine, it will be called here.
;                        Required.
;
; NOTES:
;       Requires processed GLIMPSE & BGPS cutouts generated by
;       process_glimpse_irdc.pro to be found in the
;       ./local/ubc_glimpse_proc/ directory.
;
;       Also requires the BGPS postage-stamp images produced by the
;       routine BGPS_POSTAGE.pro.
;
; MODIFICATION HISTORY:
;
;       Created:  09/08/11, TPEB -- Initial version, massaging
;                                   Erik Rosolowsky's routine
;                                   irdcdist.pro into the structure of
;                                   distance-omnibus.
;       Modified: 09/15/11, TPEB -- Added VERBOSE option, and fixed
;                                   bug in extracing LABELMAP
;                                   subsection for use.  Also, added
;                                   fix for UBC_GLIMPSE_PROC numbering
;                                   mishaps.
;       Modified: 02/27/12, TPEB -- Updated routine with all the
;                                   knowledge and work I've
;                                   done since September.
;       Modified: 03/01/12, TPEB -- Refining the prior PDF scheme, and
;                                   checking morphological matching of
;                                   objects in the IRDC modeling code
;                                   against the distance placement
;                                   done there.
;       Modified: 03/09/12, TPEB -- Added 'hyperfine' keyword to
;                                   IRDC_MORPH_FFORE for the purposes
;                                   of a plot for the paper.
;       Modified: 03/26/12, TPEB -- Modified routine to use new
;                                   fully-registered label maps
;                                   generated by
;                                   process_glimpse_irdc.pro
;       Modified: 07/02/12, TPEB -- Added constricting top-hat
;                                   function to impose on the label
;                                   map to limit spatial extent of
;                                   morphological matching region.
;       Modified: 07/04/12, TPEB -- Updated header documentation.
;       Modified: 08/03/12, TPEB -- Added ALPHA option for finding
;                                   best way to turn chisq into a
;                                   DPDF, and added memory-cleaning
;                                   lines.
;       Modified: 08/20/12, TPEB -- Moved IRDC_MORPH_FFORE to a new
;                                   routine FFORE_LBD.pro.
;       Modified: 09/05/12, TPEB -- Added IRAC BAND4 aperture
;                                   correction factor of 0.737 to
;                                   account for scattering of light
;                                   within the camera.
;       Modified: 09/06/12, TPEB -- Fixed the scattered light aperture
;                                   correction added yesterday.
;       Modified: 09/27/12, TPEB -- Shifted FFORE calculation from the
;                                   analytical calculation to using
;                                   the FFORE data cube (l,b,d).
;       Modified: 10/19/12, TPEB -- Moved the IRAC Band 4 scattering
;                                   term correction to
;                                   process_glimpse_irdo.pro so that
;                                   the postage-stamp images are
;                                   correct for all purposes.
;       Modified: 11/02/12, TPEB -- Added forcing ffore <= 1.0.
;       Modified: 11/26/12, TPEB -- Changed default alpha to 3.
;       Modified: 11/27/12, TPEB -- Changes to PostScript plots, plus
;                                   going back to alpha = 2.
;       Modified: 12/13/12, TPEB -- Adjusted the location of the
;                                   colorbars to make them line up
;                                   with the plot windows.
;       Modified: 01/14/13, TPEB -- Adding functionality for alternate
;                                   DPDF formulation, suggested by
;                                   Erik.  Also, code cleanup to
;                                   remove old testing sections.
;       Modified: 01/17/13, TPEB -- Added optional output for
;                                   alternate DPDF formulation (PROB).
;       Modified: 02/07/13, TPEB -- Made EPS plots paper-friendly.
;
;-


;; ;;===================================================================
;; ;; Subroutine for use with NEWTON & fmeas
;; FUNCTION IRDC_MORPH_DFMEAS, d, FMEAS=fmeas, L=l, B=b
  
;;   ffore = FFORE_LBD(l, b, d, /NONNATIVE, /HYPERFINE, /FROMFIT)
  
;;   RETURN,abs(fmeas - ffore)*1.d5
;; END


;;================================================
;;
;; Main Routine
;;
FUNCTION IRDC_MORPH, s, dvec, MAKE_PS=make_ps, VERBOSE=verbose, CLEM=clem, $
                     TOPHAT=tophat, TD=Td, FMEAS=fmeas, MDIST=mdist, $
                     SUFF=suff, BEAMMASK=beammask, ALPHA=alpha, $
                     UPSILON=upsilon, SILENT=silent, $
                     PDFDC=pdfdc, DIFF=diff, DOF=dof
  
  
  ;;===============================
  ;; Set clip to 1 for clipping, 0 for removing of bright GLIMPSE
  ;; pixels from within the BGPS contour for morphological matching.
  ;;
  ;; SET TO ZERO
  clip = 0b
  ;;
  ;;===============================
  
  COMMON POSTAGE_BLOCK, common_img, common_hdr
  COMMON CUBE_BLOCK, ff_cube, ff_astr, ff_dist
  
  ;; If the FFORE data cube (l,b,d) is not loaded into the COMMON
  ;;   block, do so now...
  IF total(size(ff_cube)) EQ 0 THEN omni_load_ffore
  
  
  start_t = systime(1)
  ;; Parse keywords
  verbose = KEYWORD_SET( verbose )
  silent  = KEYWORD_SET(silent)
  IF verbose && ~silent THEN $
     message,'Running IRDC_MORPH on BGPS #'+string(s.cnum,format="(I4)"),/inf
  IF n_elements(tophat) EQ 0 THEN tophat = 0
  dofdist = (n_elements(fmeas) EQ 0 || fmeas EQ 2) ? 0b : 1b
  domdist = (n_elements(mdist) EQ 0 || mdist EQ -1) ? 0b : 1b
  IF n_elements(suff) EQ 0 THEN suff=''
  beammask = KEYWORD_SET(beammask)
  IF n_elements(alpha) EQ 0 THEN alpha = 2.
  IF n_elements(Td) EQ 0 THEN Td = 20.d
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  d = (n_elements(dvec) NE 0) ? dvec : $
      dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
  
  ;; Set directory names and whatnot
  UBC_DIR = './local/ubc_glimpse_proc/'
  scnum = string(s.cnum,format="(I04)")
  
  ;; Define arrays
  prob = dblarr(n_elements(d))
  diff = dblarr(n_elements(d))
  
  ;; Check that this object has the UBC GLIMPSE PROC images!!1!
  IF ~ FILE_TEST(UBC_DIR+'bgps'+scnum+'.fits') THEN BEGIN
     message,'Source not contained in UBC_GLIMPSE_PROC...',/inf
     return,prob
  ENDIF
  
  ;; Read in the proper files for this object
  IF n_elements(common_img) EQ 0 THEN BEGIN
     bgps = readfits(UBC_DIR+'bgps'+scnum+suff+'.fits',bgpshdr, $
                     /silent) * 1.5d 
     Imir = readfits(UBC_DIR+'Imir'+scnum+suff+'.fits',imirhdr, /silent)
     irac = readfits(UBC_DIR+'smooth33arc'+scnum+suff+'.fits',irachdr, /silent)
     labl = readfits(UBC_DIR+'label'+scnum+suff+'.fits',labhdr, /silent)
  ENDIF ELSE BEGIN  ;; READ IN FROM COMMON BLOCK
     ;; message,'Running from COMMON...',/inf
     sind = WHERE(common_img.cnum EQ s.cnum, ns)
     IF ns NE 1 THEN message,'500 CCs of vodka, STAT!'
     bgps = common_img[sind].bgps
     Imir = common_img[sind].imir
     irac = common_img[sind].irac
     labl = common_img[sind].labl
     bgpshdr = common_hdr[sind].bgps
     imirhdr = common_hdr[sind].imir
     irachdr = common_hdr[sind].irac
     labhdr  = common_hdr[sind].labl
  ENDELSE
  label = labl                  ; To keep track for pixel removal
  
  extast, bgpshdr, astr
  xy2ad, 16, 16, astr, l, b     ; Center of UBC postage
  
  time1 = systime(1) - start_t
  ;;========================================================
  ;; Determine hits pixels
  labl *= (irac LE Imir)
  
  hits = where(labl, pixct)
  
  ;; If pixct=0, then no overlap, return null
  IF pixct EQ 0 THEN RETURN,prob  
  
  ;; Degrees of freedom, given that pixels are not independent.
  dof = double(pixct) / 23.8d
  
  ;;=================================================================
  ;; Calculate the optical depth map based on BGPS & Td
  IF n_elements( upsilon ) EQ 0 THEN upsilon = flux2tau_bgps( Td )
  
  taumap = upsilon * (bgps > 0.d)
  
  ;; Get statistics and the ffore from the model along the LoS
  bgpsrms = mad(bgps)
  modelrms = Imir * exp(-taumap) * upsilon * bgpsrms * (1.d)
  
  ;; Work the FFORE MAGIC!
  dfvec = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  ad2xy,s.glon_peak,s.glat_peak,ff_astr,xi,yi
  ff = (ff_cube[xi,yi,*]) < 1.0
  ffore = interpol(ff,ff_dist,dfvec) < 1.0
  
  ;; IF dofdist THEN BEGIN
  ;;    print,fmeas
  ;;    quiet = 1b
  ;;    first=1b
  ;;    functargs = {L:s.glon_peak, B:s.glat_peak, FMEAS:fmeas}
  ;;    guess = (fmeas LE 0) ? 1.d3 : !MW.R0*cos(l*!dtor)
  ;;    parinfo = {value:guess,limited:[1b,1b],limits:[0,20.d3],$
  ;;               step:100}
  ;;    dfmeas = TNMIN('IRDC_MORPH_DFMEAS',guess,STATUS=status,$
  ;;                   FUNCTARGS=functargs,/AUTODERIVATIVE,PARINFO=parinfo,$
  ;;                   ERRMSG=errmsg, QUIET=quiet, MAXITER=1000)
  ;;    IF verbose && status LE 0 && status NE -18 THEN message,$
  ;;       string(status,format="(I0,'  ')")+errmsg,/cont 
  ;;    print,'DFMEAS [kpc]: ',dfmeas/1.d3
  ;; ENDIF
  
  ;; Loop over distances!
  FOR j=0L, n_elements(d)-1 DO BEGIN
     
     ;; Model 8-um emission
     undefine,model
     model = Imir * ( ffore[j] + (1.d - ffore[j]) * exp(-taumap) )
     
     ;; \chi^2_j
     ;; Divide by number of pixels per beam!!!
     diff[j] = TOTAL(((irac-model)[hits])^2 / modelrms[hits]^2 / 23.8d)
     
  ENDFOR
  
  ;; Calculate statistics (namely degrees of freedom given pixels are
  ;;   not independent (23.80 pixels per beam)
  ;; diff /= dof
  min_diff = min(diff,j)
  
  ;; New probability
  invx = diff^(-float(alpha))
  invx /= MAX(invx)
  dof = (floor(n_elements(hits)/23.8)-1) > 1
  pdfdc = 1-chisqr_pdf(diff-min(diff),dof)
  invdx = (diff-min_diff)^(-float(alpha))
  invdx /= MAX(invdx)

  ;; yfit = mpfitpeak(d,invx,A,nterms=4)
  ;; IF ~silent THEN print,'Gauss Fit: ',A
  
  time2 = systime(1) - start_t
  ;;================================================================
  ;; If selected, generate pretty pictures
  IF KEYWORD_SET( make_ps ) THEN BEGIN
     
     COMMON VEL_BLOCK, velocity
     IF ~ exist(velocity) THEN BEGIN
        message,'Warning: COMMON VEL_BLOCK required '+$
                'for KDIST overlay',/cont
        kprob = (invx*0.d + 1.d)
        do_overlay = 0b
     ENDIF ELSE BEGIN
        this  = WHERE(velocity.cnum EQ s.cnum, nind)
        kprob = kdist_spectrum(s.glon_peak, s.glat_peak, velocity[this].v_std, $
                               velocity[this].densegas, DVEC=dvec, $
                               ERRVLSR=errvlsr, CONSTRAIN=constrain, CLEM=clem)
        kprob /= MAX(kprob)
        do_overlay = 1b
     ENDELSE
     
     ;; The model for the best-fit distance
     model = Imir * ( ffore[j] + (1.d - ffore[j]) * exp(-taumap) )
     
     ;; Clip irac image intensities for plotting purposes
     thresh = ( MEDIAN(Imir) + 1.d*MAD(irac) )
     irac = irac < thresh
     pr = set_plot_range(irac)
     
     ;; Set plot environment
     myps,'./local/output/morph'+scnum+'.eps',xsize=6.5,ysize=5.65,ct=0,/cmyk
     charsize = 0.7
     chi = cgSymbol('chi')
     sun = cgSymbol('sun')
     legcolor1 = 'WT3'
     legcolor2 = 'Black'
     
     multiplot_xm,[2,2],/xmargin,ygap=0.02,xgap=0.042,mpcharsize=charsize*0.7
     
     cgText,/norm,align=0.5,0.5,0.97,'BGPS #'+string(s.cnum,format="(I4)"),$
            charsize=charsize*1.2
     
     
     ;;=======================================
     ;; Panel 1 -- Model 8-um plot
     plotmap,model,irachdr,charsize=charsize,range=pr,XC=xc,YC=yc,ct=0
     cgContour,labl,xc,yc,levels=[0.5],thick=5,/overplot,$
               color='Crimson',label=0
     al_legend,/top,/right,['(a)'],textcolor=legcolor2,charsize=charsize*1.2,$
               box=0
     al_legend,/bottom,/left,['Synthetic'],textcolor=legcolor2,$
               charsize=charsize*1.1,box=0
     
     multiplot,/doyaxis,/doxaxis
     
     ;;=======================================
     ;; Panel 2 -- BGPS image
     plotmap,bgps>0.d,bgpshdr,charsize=charsize,ct=3,range=bpr,$
             axcolor='WT3'
     cgContour,labl,xc,yc,levels=[0.5],thick=5,/overplot,$
               color='light cyan',label=0
     al_legend,/top,/right,['(b)'],textcolor=legcolor1,charsize=charsize*1.2,$
               box=0
     al_legend,/bottom,/left,['BGPS'],textcolor=legcolor1,$
               charsize=charsize*1.1,box=0
     bywin = !y.window
     bxwin = !x.window
     
     multiplot,/doxaxis,/doyaxis
    
     ;;=======================================
     ;; Panel 3 -- Smoothed GLIMPSE image
     plotmap,irac,irachdr,charsize=charsize,range=pr,ct=0
     cgContour,labl,xc,yc,levels=[0.5],thick=5,/overplot,$
               color='Crimson',label=0
     al_legend,/top,/right,['(c)'],textcolor=legcolor2,charsize=charsize*1.2,$
               box=0
     gywin = !y.window
     al_legend,/bottom,/left,['GLIMPSE'],textcolor=legcolor2,$
               charsize=charsize*1.1,box=0
     
     multiplot,/doxaxis,/doyaxis
     
     ;;=======================================
     ;; Panel 4 -- CHISQ plot
     cgPlot, d/1.d3, invx, xtit='Heliocentric Distance  [kpc]',yr=[0,1.05],$
             /yst,/nodata,ytit='Relative Probability',ytickformat="(F0.1)",$
             charsize=charsize,position=[bxwin[0],gywin[0],bxwin[1],gywin[1]]
     
     ;; Plot D_tan & ffore(d)
     tandist = !MW.R0*cos(l * !dtor)/1.d3
     cgOplot,d/1.d3,ffore,color='BLK6',linestyle=1,thick=3
     vline,tandist,color='Green',thick=6
     cgText,tandist+0.2,0.98,['d!dtan!n'],charsize=charsize,color='Green'
     
     IF dofdist THEN BEGIN
        vline,dfmeas/1.d3,color='BLU6',thick=4
        cgText,12,1.2,charsize=charsize,'f!dfore,meas!n = '+$
               string(fmeas,format="(F0.3)"),color='BLU6'
     ENDIF     
     IF domdist THEN vline,mdist,color='RYB3',thick=4
     
     ;; Molecular gas prior
     h2_prob = prob_lonlat(s)
     h2_prob /= max(h2_prob)
     cgOplot,d/1.d3,h2_prob,thick=3,color='BLU5',linestyle=3
     
     ;; Plot Morph Match DPDF
     cgOplot,d/1.d3, invx,thick=6
     
     ;; cgOplot,d/1.d3, pdfdc, thick=6, color='dodger blue'
     ;; cgOplot,d/1.d3, invdx, thick=6, color='orchid'
     
     IF do_overlay THEN BEGIN
        ;; oplot,d/1.d3, kprob, thick=6, color=probcolor(3)
        totprob = (kprob * invx * h2_prob) / MAX(kprob * invx * h2_prob)
        cgOplot,d/1.d3, totprob, thick=6, color='Crimson'
        ;; tp2 = (kprob * pdfdc * h2_prob) / MAX(kprob * pdfdc * h2_prob)
        ;; cgOplot,d/1.d3, tp2, thick=6, color='Brown'
        
        ;; ind = WHERE(totprob EQ max(totprob), nind)
        ;; IF nind NE 1 THEN tag = 'A' ELSE BEGIN
        ;;    IF d[ind]/1.d3 LE tandist THEN tag = 'N' ELSE tag = 'F'
        ;; ENDELSE
        ;; al_legend,tag,pos=[16,0.2],textcolor='PUR7',box=0
        
     ENDIF
     
     ;; fi = where(d/1.d3 GE tandist, nf)
     ;; ni = where(d/1.d3 LT tandist, nn)
     ;; max1 = max(totprob,j1)
     ;; max2 = max(tp2,j2)
     ;; totprob /= total(totprob)
     ;; tp2     /= total(tp2)
     ;; pch1 = d[j1]/1.d3 GE tandist ? total(totprob[fi]) : total(totprob[ni])
     ;; pch2 = d[j2]/1.d3 GE tandist ? total(tp2[fi]) : total(tp2[ni])
     ;; cgText,14,0.8,['P!dML!n = '+string(pch1,format="(F0.3)")],$
     ;;        color=probcolor(2),charsize=charsize
     ;; cgText,14,0.7,['P!dML!n = '+string(pch2,format="(F0.3)")],$
     ;;        color='Brown',charsize=charsize
     ;; ;; cgText,14,0.6,['NF: '+string(nf,format="(I0)")],color='green',$
     ;; ;;        charsize=charsize
     ;; ;; cgText,14,0.5,['NN: '+string(nn,format="(I0)")],color='green',$
     ;; ;;        charsize=charsize
     ;; cgText,14,0.6,[cgSymbol('chi')+'!u2!dmin!n = '+$
     ;;                string(min_diff,format="(F0.1)")],color='green',$
     ;;        charsize=charsize
     ;; cgText,14,0.5,['dof = '+string(dof,format="(I0)")],color='green',$
     ;;        charsize=charsize
     ;; cgText,14,0.4,['d!dML!n = '+string(d[j1]/1.d3,format="(F0.2)")+' kpc'],$
     ;;        color=probcolor(2),charsize=charsize*0.8
     ;; cgText,14,0.3,['d!dML!n = '+string(d[j2]/1.d3,format="(F0.2)")+' kpc'],$
     ;;        color='Brown',charsize=charsize*0.8

    
     axis,xaxis=0,xthick=3,/xst,xtickformat='blank_axis'
     al_legend,/top,/right,['(d)'],textcolor=legcolor2,charsize=charsize*1.2,$
               box=0
     
     cgLoadct,0,/silent
     cgcolorbar,range=pr,position=[0.905,gywin[0],0.925,gywin[1]],$
                divisions=0,title='GLIMPSE Intensity  [MJy sr!u-1!n]',$
                charsize=charsize,/vertical,/right
     
     cgLoadct,3,/silent
     cgcolorbar,range=bpr,position=[0.905,bywin[0],0.925,bywin[1]],$
                divisions=0,title='BGPS Flux Density  [Jy beam!u-1!n]',$
                charsize=charsize,/vertical,/right
     myps,/done,/mp

  ENDIF  ;; End of postscript-creating IF block
  
  ;; Clean Memory
  undefine,bgps,Imir,irac,label,labl
  
  time3 = systime(1) - start_t
  ;; print,time1,time2,time3
  
  junk = CHECK_MATH()
  ;; RETURN, pdfdc ;; Return the inverse chi-squared function as the prior PDF
  RETURN,invx  ;; Return the inverse chi-squared function as the prior PDF
END
